Tarea3

Minimos cuadrados

Plantear el problema de minimos cuadrados como un problema de minimizacion:

Tenemos que \(Y=X\beta + \epsilon\), por lo que donde \(\epsilon\) es el error, por lo que podemos plantear el problema como la minimizaci´ón de este error \[min_{ } ||\epsilon||^2 \] que es equivalente a \[min_{\beta} ||Y-X\beta||^2\] Para encontrar el estimador \(\hat{\beta}\) es necesario derivar e igualar a 0: Derivando respecto a \(\beta\): \[-2X^t (Y-X\hat{\beta})=0\] \[2X^tY=2X^tX\hat{\beta}\] \[\hat{\beta}=(X^tX)^{-1}X^tY\] Es el \(\hat{\beta}\) que soluciona el problema de minimización.

Con esto podemos ver que para cada: \[\hat{\beta}_{i}=\frac{\sum_{i}x_{i}y_{i}}{\sum_{i}x_{i}^2}\] La cuál es una composición lineal de los parametros. Desde el planteamiento estamos buscando que el resultado sea lineal en los parametros, no en los datos por lo que este metodo serviria para poder aproximar polinomios de la orma \(Y=X^2\)

¿En qué se relaciona esto con un problema de proyección?

Recordemos que al encontrar la proyección de un vector \(Proy_{X}(y)= \frac{<x,y>}{||x||^2}y\) el error de proyeccion \(Y- Proy_{x}(Y)\) será ortogonal a la proyección. En el caso del problema de minimizar el error, nosotros estamos proyectando la variable dependiente Y sobre el espacio formado por nuestros datos X. y buscamos las combinaciones sobre X que hagan este error mas pequeño.

Viendolo como el teorema de Pitagoras el cu´ál nos dice que el cuadrado de la hipotenisa de un triangulo es igual a la suma de los otros dos lados \(h^2=x^2+y^2\) si el angulo formado entre \(x\) y \(y\) es de 90 grados, en este caso \(H=Y\) las variables dependientes, \(x=Proy_{x}(Y)\) y \(y= error\), entonces la manera en la que se puede cumplir el teorema es que el error sea ortogonal a la proyeccion (angulo de 90 grados) sino el error seráa mas grande.

¿Que logramos al agregar una columna de unos en la matriz X?

Al definir una columna de unos estamos cambiando un poco el problema. Antes de agregar la columna: \[y_{j}=\beta_{1}x_{1j}+...+\beta_{n}x_{nj}\] Al definir la columa de unos estamos dando cierta holgura al modelo ya que al agregar un parametro \(\alpha\) no estamos forzando a que la aproximacion pase por el origen. \[y_{j}=\alpha+\beta_{1}x_{1j}+...+\beta_{n}x_{nj}\]

Plantear el problema de regresion como un problema de estadistica con errores….

¿Cual es la funcion de verosimilitud del problema anterior?

Planteando el problema como \[Y= X\beta +\epsilon\] con \[\epsilon \sim N(0,\sigma^2I_{n})\]

podemos probar que \[Y \sim N(X\beta,\sigma^2I_{n})\]

Demostracion: \[Y= X\beta +\epsilon\] \[\epsilon=Y-X\beta\]

1. Media: \[E(\epsilon)=E(Y)-E(X\beta)=0\] \[E(Y)=E(X\beta)\] \[E(Y)=X\beta\]
2. Varianza: \[Var(\epsilon)=Var(Y)+Var(X\beta)-2Cov(Y,X\beta)=\sigma^2I_{n}\] donde: \[Var(X\beta)=0\] y \[Cov(Y,X\beta)=E((Y-E(Y)(X\beta-E(X\beta)))\] \[=E((Y-X\beta)(X\beta-X\beta))=0\] entonces: \[Var(\epsilon)=Var(Y)=\sigma^2I_{n}\]

Ahora, falta demostrar que Y es Normal: Por propiedades de la Normal como \(\epsilon \sim N(0,\sigma^2I_{n})\) \[\epsilon + X\beta \sim N(X\beta, \sigma^2I_{n})\] (Esto era m´ás rapido pero bueno…)

Entonces con esto podemos escribir la funcion de verosimilitud de Y como: \[ L(\beta, \sigma^2)=\prod \frac{1}{\sqrt{2\pi\sigma^2}} exp(-\frac{(y-X\beta)^2}{2\sigma^2})\]

\[=(\frac {1}{\sqrt{2\pi\sigma^2}})^n exp (- \frac{1}{2\sigma} (y-X\beta)^t(y-X\beta))\]

A la cual le sacamos Logaritmo para despu´és poder buscar los parametros que maximicen la funcion: \[Log(L(\beta, \sigma^2))=-\frac{n}{2}ln(2\pi)-\frac{n}{2}ln(\sigma^2)-\frac{1}{2\sigma^2}(y-X\beta)^t(y-X\beta)\] Ahora hay que maximizar esa funcion sobre \(\beta\) y \(\sigma^2\) Respecto de \(\beta\) \[\frac{dL}{d\beta}=\frac{1}{\sigma^2}(y-X\beta)^tX=0\] \[Xy=X^tX\beta\] \[\hat{\beta}=(X^tX)^{-1}Xy\]

Respecto de \(\sigma^2\) \[\frac{dL}{d\sigma^2}= - \frac{n}{2\sigma^2}+\frac{1}{2\sigma^4}(y-X\beta)^t(y-X\beta)=0\] \[\sigma^2=\frac{1}{n} (y-X\beta)^t(y-X\beta)\] que puede resolverse usando \(\hat{\beta}\) que encontramos en la solucion anterior.

Como pudimos ver el resultado de maxima verosimilitud \[\hat{\beta}=(X^tX)^{-1}Xy\] Es igual al resultado de \(\hat{\beta}\) de minimos cuadrados

El teorema de Gauss Markov

El teorema de Gauss Markov plantea que un modelo de regresion lineal en el que se cumple que:

1. \(E(\epsilon)=0\)
2. Los errores no estan correlacionados entre si
3. \(Var(\epsilon_{i}=\sigma^2\) para todo i, los errores tienen la misma varianza

Entonces el mejor estimador lineal insengado de los oeficientes \(\beta\) es el resultante del problema de minimizar los errores al cuadrado. si existe.

Parte Aplicada

summary(reg_precio)

Call:
lm(formula = precio ~ carat + depth + table + x + y + z, data = diamonds_num)

Residuals:
     Min       1Q   Median       3Q      Max 
-23878.2   -615.0    -50.7    347.9  12759.2 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 20849.316    447.562  46.584  < 2e-16 ***
carat       10686.309     63.201 169.085  < 2e-16 ***
depth        -203.154      5.504 -36.910  < 2e-16 ***
table        -102.446      3.084 -33.216  < 2e-16 ***
x           -1315.668     43.070 -30.547  < 2e-16 ***
y              66.322     25.523   2.599  0.00937 ** 
z              41.628     44.305   0.940  0.34744    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1497 on 53933 degrees of freedom
Multiple R-squared:  0.8592,    Adjusted R-squared:  0.8592 
F-statistic: 5.486e+04 on 6 and 53933 DF,  p-value: < 2.2e-16

Que tan bien ajusta el modelo depende de la correlacion que tenga cada una de las variables explicativas sobre el precio, esto lo podemos ver a trav´és de:

cor(diamonds_num)

           price      carat       depth      table           x           y          z
price  1.0000000 0.92159130 -0.01064740  0.1271339  0.88443516  0.86542090 0.86124944
carat  0.9215913 1.00000000  0.02822431  0.1816175  0.97509423  0.95172220 0.95338738
depth -0.0106474 0.02822431  1.00000000 -0.2957785 -0.02528925 -0.02934067 0.09492388
table  0.1271339 0.18161755 -0.29577852  1.0000000  0.19534428  0.18376015 0.15092869
x      0.8844352 0.97509423 -0.02528925  0.1953443  1.00000000  0.97470148 0.97077180
y      0.8654209 0.95172220 -0.02934067  0.1837601  0.97470148  1.00000000 0.95200572
z      0.8612494 0.95338738  0.09492388  0.1509287  0.97077180  0.95200572 1.00000000

Como podemos ver en la tabla de correlaciones, precio est´á correlacionada altamente con “carat”, “x”,“y”,“z” y tienen cierta correlacion con death y table, esto me dice que las variables escogidas sirven para explicar el precio. Podemos ver esta relacion de manera grafica a continuacion en el siguiente Diagrama de dispersion :

library(GGally)
library(ggplot2)
my_fn <- function(data, mapping, ...){
  p <- ggplot(data = data, mapping = mapping) + 
    geom_point() + 
    geom_smooth(method=loess, fill="red", color="red", ...) +
    geom_smooth(method=lm, fill="blue", color="blue", ...)
  p
}
g = ggpairs(diamonds_num,columns = 1:7, lower = list(continuous = my_fn))
g

 plot: [1,1] [==-------------------------------------------------------------------------------------------------------]  2% est: 0s 
 plot: [1,2] [====-----------------------------------------------------------------------------------------------------]  4% est: 4s 
 plot: [1,3] [======---------------------------------------------------------------------------------------------------]  6% est: 4s 
 plot: [1,4] [=========------------------------------------------------------------------------------------------------]  8% est: 4s 
 plot: [1,5] [===========----------------------------------------------------------------------------------------------] 10% est: 4s 
 plot: [1,6] [=============--------------------------------------------------------------------------------------------] 12% est: 4s 
 plot: [1,7] [===============------------------------------------------------------------------------------------------] 14% est: 4s 
 plot: [2,1] [=================----------------------------------------------------------------------------------------] 16% est: 4s 
 plot: [2,2] [===================--------------------------------------------------------------------------------------] 18% est: 2m 
 plot: [2,3] [=====================------------------------------------------------------------------------------------] 20% est: 2m 
 plot: [2,4] [========================---------------------------------------------------------------------------------] 22% est: 2m 
 plot: [2,5] [==========================-------------------------------------------------------------------------------] 24% est: 1m 
 plot: [2,6] [============================-----------------------------------------------------------------------------] 27% est: 1m 
 plot: [2,7] [==============================---------------------------------------------------------------------------] 29% est: 1m 
 plot: [3,1] [================================-------------------------------------------------------------------------] 31% est: 1m 
 plot: [3,2] [==================================-----------------------------------------------------------------------] 33% est: 2m 
 plot: [3,3] [====================================---------------------------------------------------------------------] 35% est: 3m 
 plot: [3,4] [=======================================------------------------------------------------------------------] 37% est: 2m 
 plot: [3,5] [=========================================----------------------------------------------------------------] 39% est: 2m 
 plot: [3,6] [===========================================--------------------------------------------------------------] 41% est: 2m 
 plot: [3,7] [=============================================------------------------------------------------------------] 43% est: 2m 
 plot: [4,1] [===============================================----------------------------------------------------------] 45% est: 2m 
 plot: [4,2] [=================================================--------------------------------------------------------] 47% est: 2m 
 plot: [4,3] [===================================================------------------------------------------------------] 49% est: 2m 
 plot: [4,4] [======================================================---------------------------------------------------] 51% est: 3m 
 plot: [4,5] [========================================================-------------------------------------------------] 53% est: 2m 
 plot: [4,6] [==========================================================-----------------------------------------------] 55% est: 2m 
 plot: [4,7] [============================================================---------------------------------------------] 57% est: 2m 
 plot: [5,1] [==============================================================-------------------------------------------] 59% est: 2m 
 plot: [5,2] [================================================================-----------------------------------------] 61% est: 2m 
 plot: [5,3] [==================================================================---------------------------------------] 63% est: 2m 
 plot: [5,4] [=====================================================================------------------------------------] 65% est: 2m 
 plot: [5,5] [=======================================================================----------------------------------] 67% est: 2m 
 plot: [5,6] [=========================================================================--------------------------------] 69% est: 2m 
 plot: [5,7] [===========================================================================------------------------------] 71% est: 2m 
 plot: [6,1] [=============================================================================----------------------------] 73% est: 2m 
 plot: [6,2] [===============================================================================--------------------------] 76% est: 2m 
 plot: [6,3] [=================================================================================------------------------] 78% est: 2m 
 plot: [6,4] [====================================================================================---------------------] 80% est: 2m 
 plot: [6,5] [======================================================================================-------------------] 82% est: 1m 
 plot: [6,6] [========================================================================================-----------------] 84% est: 1m 
 plot: [6,7] [==========================================================================================---------------] 86% est: 1m 
 plot: [7,1] [============================================================================================-------------] 88% est: 1m 
 plot: [7,2] [==============================================================================================-----------] 90% est:50s 
 plot: [7,3] [================================================================================================---------] 92% est:41s 
 plot: [7,4] [===================================================================================================------] 94% est:32s 
 plot: [7,5] [=====================================================================================================----] 96% est:22s 
 plot: [7,6] [=======================================================================================================--] 98% est:11s 

Graficamente así se ven los valores del modelo contra los observados

data<-as.data.frame(cbind(reg_precio$fitted.values,diamonds_num$price))
p <- plot_ly(data = data, x = ~data[,1], y = ~data[,2]) %>%
  layout(title="Valores ajustadosvs Valores observados")
p

No trace type specified:
  Based on info supplied, a 'scatter' trace seems appropriate.
  Read more about this trace type -> https://plot.ly/r/reference/#scatter
No scatter mode specifed:
  Setting the mode to markers
  Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode
No trace type specified:
  Based on info supplied, a 'scatter' trace seems appropriate.
  Read more about this trace type -> https://plot.ly/r/reference/#scatter
No scatter mode specifed:
  Setting the mode to markers
  Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode

Podemos ver que como el modelo tiene una \(R^2 ~ .8592\) por lo que podemos decir que el modelo esta bien especificao ya que explica el aproximadamente el 85% de la varianza con las variables que estamos usando. nuestro error estandar residual \(\sigma=1497\) por lo que \(sigma^2=1497^2\)

¿Cual es el angulo entre Y y \(\hat{Y}\)

\[arcos(R^2)=\theta\] En este caso en particular \(R^2=.8592\) por lo que

acos(.8592)

[1] 0.5370923

Haciendo una funcion de logverosimilitud

En caso de que podamos aproximar con una normal tenemos lo siguiente:

mu<-X%*%beta

Error in X %*% beta : non-conformable arguments

En nuestro caso en particular:

X_reg<-cbind(rep(1,length(y)), diamonds_num[,2:7])
y_reg<-diamonds_num[,1]
betasigma<-c(0,0,0,0,0,0,0,80) 
optim(par=betasigma,lik1, X=X_reg, y=y_reg)

Error in optim(par = betasigma, lik1, X = X_reg, y = y_reg) : 
  objective function in optim evaluates to length 0 not 1